Molecular Modeling of Indeno [1, 2-b] Quinoline-9, 11-Diones as Cytotoxic Agents

Deoxyribonucleic acid (DNA) is an important molecular target for anti-cancer agents due to its involvement in gene expression and protein synthesis which are fundamental steps in cell division and growth. A number of antineoplastic agents interfere with DNA and hence disturb the cell cycle. Compounds including planar aromatic rings are privileged scaffolds in binding to DNA. This characteristic is mainly arisen from the fact that such structural feature may be appropriate to insert between the base pairs of the DNA double helix and produce relatively stable non-covalent complexes. Besides π-π stacking interactions, binding to the DNA molecule might be intensified through H-bond interactions of heterocyclic rings. In the present contribution, a series of experimentally validated cytotoxic indeno[1,2-b]quinoline-9,11-diones (1-12) and their aromatized analogues (13-21) developed in our group were subjected to docking and molecular dynamics simulations to elucidate their most probable binding modes with DNA.


Introduction
Malignant tumors are responsible for cancer (1) and can be broken away and spread out to other tissues of the body (2). In spite of recent advances in treatment strategies, cancer is still one of the major causes of death worldwide (3). Many anti-cancer agents interact directly with DNA or prevent its processing through interaction with a definite enzyme (4, 5). Regarding to this, development of DNA targeting anti-cancer drugs is one of the considerable interests within pharmaceutical scientists.
Generally, ligand-DNA interactions can be categorized into two major classes; intercalation and groove binding (6). In intercalation mechanism, a planar aromatic/ polyaromatic ligand with sufficiently large surface area inserts into the space between the base pairs of DNA double helix. The driving forces for intercalation are primarily π-π stacking interactions, H-bonds and ionic interactions between ligand and base pairs (7). It has been postulated that intercalation reduces the DNA helical twist and elongates the DNA, hence disturbing its normal function (8).
Minor groove-binding agents usually fit to the helical turn of the DNA grooves via Van der Waals or hydrogen bond interactions (5). This structural complementary may be attributed to the torsional freedom of interfering ligands that provide compatibility with a twisted helical turn of the DNA grooves. Such molecular pattern may be envisaged as several unfused aromatic rings being connected through rotatable bonds (9). Major groove-binders have been less reported due to the requirement of larger molecular sizes in order to be fitted into the DNA. Major groove binding pattern have frequently been observed in peptides and proteins (10).
X-ray crystallographic structures of DNAligand complexes (http://www.rcsb.org) provided studies on the molecular mechanism and structure binding relationship of DNA interfering agents. Such structural data facilitated the performance of target based hit/lead design strategies toward developing new DNA binding anti-neoplastic agents.
In our previous work, a series of new imidazole substituted indeno [1,2-b] quinoline-9,11-dione derivatives were synthesized and evaluated for their cell-based cytotoxic effects on HeLa, LS180, Jurkat and MCF-7 human cancer cell lines using MTT assay (11). Prepared compounds exhibited no (IC 16 >100 µM) to good (IC 16 = 0.7 µM in Jurkat cell lines) anti-tumoral activities within evaluated cell lines. Due to the uncertainty of results in IC 50 scale, IC 16 values were used for our modeling studies. Prepared compounds had rigid aromatic/heteroaromatic rings constructed around a dihydropyridine (DHP) or pyridine core, an appropriate pattern for DNA intercalation and H-bonding. In the present contribution, docking and molecular dynamics methods were applied to elucidate the most probable DNA binding mode of cytotoxic indeno[1,2-b]quinoline-9,11-dione derivatives.

Preparation of oligonucleotide structures
Crystal structures of two DNA-ligand complexes with PDB IDs 1D32 (12) and 102D (13) were retrieved from the Protein Data Bank as representatives of macromolecular templates (14). For preparation of required DNA targets, co-crystallographic ligands DIT and TNT were respectively removed from the PDB files 1D32 and 102D and all the crystallographic water molecules were also eliminated from each PDB file. Missing hydrogens were added, non-polar hydrogens were merged to their attached carbon atoms and kollman charges were dedicated to the oligonucleotide receptors. All the pre-processing steps for receptor files were performed via AutoDock Tools program (ADT) (15,16). In addition; all the oligonucleotide structures were optimized to minimize the crystallographic induced clashes via steepest descent method by Gromacs package (17).

Preparation of ligand structures
Twenty-one DHP derivatives were categorized into indeno [1, 2-b]quinoline-9,11diones and their aromatized derivatives. 2D chemical structures were rendered by Marvin Sketch online chemical editor (18) and relevant SMILES strings were used to generate the 3D structures by free online 3D conformation generator Frog 1.0 software (19).

Molecular docking
Advanced docking package AutoDock4.2 was used to dock the screened molecules. For this purpose, 3D optimized structures of indeno,1] 2-b]quinoline-9,11-diones and their aromatized derivatives in PDB format were used as input files for the AutoDock Tools (ADT) (15, 16). All the pre-processing steps for ligand structures were performed according to the previous procedure (20). Prepared molecules were docked into the 3D structure of oligonucleotides extracted from the crystallographic files (1D32 and 102D).
Modeling of the DNA-ligand interactions were conducted within desired number of genetic algorithm runs (50 GA runs) and energy evaluations (10 6 and 5×10 5 for 1D32 and 102D, respectively) on the basis of docking validation study. Other AutoDock operating parameters were set at their default values. A map was assigned to each atom type of the ligand and DNA. Relevant maps were estimated using AutoGrid module (part of the AutoDock package). A size of the grid was set to the 50×50×50 Ǻ 3 (distributed in the x, y, and z directions) and centered on the macromolecule with a spacing of 0.375 Ǻ to include the space occupied by the DNA double helix.
LIGPLOT program was used to generate the 2D binding interaction maps between DHP structures and DNA double helix (21).

Molecular dynamics
All MD simulations were performed using GROMACS v.4.6.5 (22) applying Amber99SB (23) force field. The compatible topologies of the ligands were prepared using the General Amber Force Field (GAFF) (24) and ANTECHAMBER suite of programs (25, 26), applying the semiempirical AM1BCC charges (27,28). DNAligand complexes were prepared based on the best docking pose associated with largest binding energy.
DNA-ligand complex was centered in a dodecahedron box filled with TIP3P water (29). Enough number of Na + atoms was added to achieve electro-neutrality. A distance of 1 nm was set between the box walls and each edge of complex. For all of the molecular dynamic simulations, the Leonard-Jones and electrostatic interactions was computed using a cutoff value of 1.2 nm. The leapfrog algorithm (30) using a 1fs time-step was used as an integrator for the equation of motion. All covalent bonds were constrained using LINCS algorithm (31).
Actual simulation was set up by minimizing the systems to make sure that introduced steric clashes during the preparing process were removed. During minimization, maximum number of 5000 steps was applied using steepest descent algorithm. After initial relaxation, the system was equilibrated for 200 ps using the NVT and NPT ensemble respectively. In NVT ensemble (constant number, volume, and temperature) temperature was set at 300 K using a velocity-rescale thermostat (modified Berendsen thermostat) (32) with a coupling constant of 0.1 ps. After temperature stabled at 300 K, NPT ensemble (constant number, pressure, and temperature) was performed with an isotropic pressure of 1 bar and a coupling constant of 2 ps by using a Parrinello-Rahman barostat (33).
Position restraints were applied to the ligands and DNA during the NVT and NPT ensemble.
Long range electrostatic interactions were computed utilizing the particle mesh Ewald method (34). Equilibration was followed by a 50 ns production run using the NPT ensemble with the velocity-rescale thermostat (32) and the Parrinello-Rahman barostat (33). Evaluation of MD simulations was performed by extracting data from the trajectory files produced during the simulations.

Internal validation
The internal validation step was performed via extracting the structures of cocrystallized ligands ditercalinium (DIT) and propamidine (TNT) and re-docking them into the original crystallographic files of the d(CGCAAATTTGCG) 2 and d (CGCG) 2 oligonucleotides. Validation of the AutoDock method for each crystallographic structure was interpreted in terms of root mean square deviation (RMSD) of the Cartesian coordinates of the atoms of the ligand in the docked and crystallographic poses (RMSD < 2Ǻ). It should be noted that top-ranked AutoDock clusters in the output files were supported by high conformational population. Obtained results indicated that the parameters set for AutoDock were reliable for reproducing X-ray complex structures.
For more elucidation, lowest energy poses of cognate drugs (propamidine and ditercalinium) in docked and crystallographic states are shown in Figure 1.

Binding model selection
AutoDock is a widely used docking package in computational medicinal chemistry and it is believed to offer relatively logical outcomes within several calculations (35). To further elucidate the interactions of DNA with indeno [1,2-b]quinoline-9,11-diones (1-21), we used AutoDock4.2 software to dock all the relevant structures into the d (CGCG) 2 and d (CGCAAATTTGCG) 2 oligonucleotides of DNA double helix.
Structure elucidation of the evaluated indeno [1, 2-b]quinoline-9,11-diones showed that these compounds possessed aromatic/heteroaromatic rings constructed around a DHP/pyridine core. In our opinion, such molecular pattern might be appropriate for intercalation and H-bond interactions with DNA double helix. To explain more, possible structure DNA binding relationship of indeno [1, 2-b] quinoline-9,11diones is depicted below (Scheme 1). Observed binding potentialities persuaded us toward performing modeling studies with the aim of proposing novel DNA interfering agents.
Chemical structures of screened ligands are summarized in Figure 2. All of the molecules could be interpreted as drug-like chemical entities and for this reason they might be good candidates for further drug development strategies (36). In addition, selected DHP structures possessed 4 that top-ranked AutoDock clusters in the output files were supported by high conformational population. Obtained results indicated that the parameters set for AutoDock were reliable for reproducing X-ray complex structures.
For more elucidation, lowest energy poses of cognate drugs (propamidine and ditercalinium) in docked and crystallographic states are shown in Figure 1. Structure elucidation of the evaluated indeno [1, 2-b]quinoline-9,11-diones showed that these compounds possessed aromatic/heteroaromatic rings constructed around a DHP/pyridine core. In our opinion, such molecular pattern might be appropriate for intercalation and H-bond interactions with DNA double helix. To explain more, possible structure DNA binding relationship of indeno [1, 2-b] quinoline-9,11-diones is depicted below (Scheme 1).
Observed binding potentialities persuaded us toward performing modeling studies with the aim of proposing novel DNA interfering agents. Scheme 1. Possible relationship between chemical structure and DNA binding ability of evaluated indeno [1, 2-b] quinoline-9, 11-diones, as can be predicted from the model, both hydrophobic and hydrogen interactions might be expected on the basis of represented scaffold template.
Chemical structures of screened ligands are summarized in Figure 2. All of the molecules could be interpreted as drug-like chemical entities and for this reason they might be good candidates for further drug development strategies (36). In addition, selected DHP structures possessed low number of active torsions (rotatable bonds) which proposed them as suitable candidates for docking simulation due to lower free torsional energies. Scheme 1. Possible relationship between chemical structure and DNA binding ability of evaluated indeno [1, 2-b] quinoline-9, 11-diones, as can be predicted from the model, both hydrophobic and hydrogen interactions might be expected on the basis of represented scaffold template. low number of active torsions (rotatable bonds) which proposed them as suitable candidates for docking simulation due to lower free torsional energies.

Binding affinity prediction
Predicted docked energies can be reported in terms of final docked energy (FDE) which is the sum of final intermolecular energy and final internal energy of ligand or estimated free energy change of binding (ΔG) that is the sum of final intermolecular energy and torsional free energy penalty (37).

Binding affinity prediction
Predicted docked energies can be reported in terms of final docked energy (FDE) which is the sum of final intermolecular energy and final internal energy of ligand or estimated free energy change of binding (ΔG) that is the sum of final intermolecular energy and torsional free energy penalty (37). 1D32 and 102D macromolecular templates, respectively. For more clarification, AutoDock binding affinities together with probability factors in terms of number of clusters and top-ranked cluster population in AutoDock output files are summarized in Table 1.
Comparative analysis of experimental cytotoxicity values and calculated binding affinities might be interpreted via following rationales: (1) Nearly all aromatized indeno [1, 2-b] quinoline-9,11-diones were associated with higher binding energies when compared to their non-aromatized analogues ( Figure  3). One of the possible explanations for such observation might be attributed to the superior π-π stacking interactions with oligonucleotides of DNA double helix in aromatized compounds. Such observation could not be confirmed by experimental cytotoxic activities in Jurkat cell lines (Such rationalization may not be offered for other cell lines due to the uncertainty of cytotoxic data). This trend might be explained by this point that the cytotoxic effect of aromatized scaffolds in Jurkat cell lines is not related to the DNA intercalation.
(3) FDE values over -10 kcal/mol were contributed to compounds 3, 9, 14, 16 & 20 in 102D model. Experimental results showed that except for compound 14 which was inactive in HeLa, LS180 and MCF7 cell lines, the other four compounds showed good anti-tumoral effects within tested cells. In relation to this, modeling studies confirmed that the incorporation of electron-withdrawing substituents such as chlorine atoms to the imidazole moiety of indeno [1, 2-b] quinoline-9, 11-dions might be associated with higher cytotoxic effects.
(4) Linear regression analysis of calculated binding affinities and experimental cytotoxicity values showed that maximum correlation could be achieved in the case of HeLa cells (r0.74 = ; r 2 0.56 = for experimental IC 16 versus predicted K i values in 102D system) ( Figure 4). No correlation could be found between predicted DNA affinities and cytotoxic activities in Jurkat cell lines.

Ligand-DNA binding interactions
With the aim of further elucidating ligand-DNA binding interactions, LIGPLOT program was used to generate the binding interactions between DHPs and DNA. Selective hydrophobic interactions were depicted using the hydrophobic bond module of Ligplot (Supplementary material 1).
(3) FDE values over -10 kcal/mol were contributed to compounds 3, 9, 14, 16 & 20 in 102D model. Experimental results showed that except for compound 14 which was inactive in HeLa, LS180 and MCF7 cell lines, the other four compounds showed good anti-tumoral effects within tested cells. In relation to this, modeling studies confirmed that the incorporation of electron-withdrawing substituents such as chlorine atoms to the imidazole moiety of indeno [1, 2-b] quinoline-9, 11-dions might be associated with higher cytotoxic effects.
(4) Linear regression analysis of calculated binding affinities and experimental cytotoxicity values showed that maximum correlation could be achieved in the case of HeLa cells (r = 0.74; r 2 = 0.56 for experimental IC 16 versus predicted K i values in 102D system) ( Figure 4). No correlation could be found between predicted DNA affinities and cytotoxic activities in Jurkat cell lines. 2) Model 102D provided more hydrophobic participants when compared to 1D32. This might be expected since the residual constituents of 102D template included additional adenine and thymine bases with thymine being the most hydrophobic interacted residue.
3) All the compounds exhibited hydrophobic contacts with a number of DNA nucleobases in  7) Unlike the model 1D32, in model 102D all the non-aromatized compounds participated in key H-bond interactions via DHP NH site. Compounds 6-10 made H-bonds with sugar ring oxygen of nucleotide while binding template for compounds 1-3 was through thymine carbonyl oxygen. 8) N-phenyl substituted compounds (4, 5, 11 and 12) contributed to additional hydrophobic interactions with DC7 and DG8 nucleobases of DNA (1D32). Such hydrophobic contacts were observed within DA18, DT8, DT9 and DT19 residues in model 102D.

Molecular dynamic simulations
Dynamic stabilities of two indeno [1, 2-b] quinoline-9,11-dione/DNA complexes were evaluated using MD simulations. Considering probable binding mechanisms, i.e. intercalation (1D32) and minor groove binding (102D), the ligand with highest binding energy was selected for further simulations in each case.
Due to the planar structure of Compound 16, it was expected that this compound could bind to the DNA structure through intercalating between DG and DC nucleobases. Moreover; presence of two methyl groups at position 7 of quinoline moiety opposed steric hindrance on compound 18 and made its interaction more favorable within DNA minor groove when compared to the intercalation between DNA nucleobases.
After solvation of the DNA-ligand systems in a dodecahedron box with TIP3P water and adding enough number of Na + ions to achieve electro-neutrality, 5000 minimization steps was performed using the steepest descent method. Soaking of DNA-ligand system was carried out by 100 ps constrained NVT followed by 100 ps NPT ensembles. Following system equilibration, 50 ns MD production was performed in each case without any constraint.
Total energy, temperature and root mean square deviation (RMSD) were assessed to confirm the stability of trajectories. The average temperature during MD simulations at 300 K was found to be 299.0 K (± 3.6) and 298.2 K (± 2.3) in the case of models 1D32 and 102D, respectively. The temperatures of the both systems were stable for the whole simulation time. Fluctuation of temperature and energy is depicted in Figure 5. Obtained results showed that the conservation of energy over 50 ns MD simulations was convinced in both systems.
The stability of the MD simulation was determined in terms of departures and fluctuations from the initial ligand-DNA structure. The time evolution of the ligand (all atoms) RMSDs was recorded as a function of time. The RMSD of two ligands with regard to initial conformations was illustrated in Figure 6.
As can be seen in Figure 7, both ligands had overall stability in binding state. Compound 16 bound to DNA through intercalation between DC and DG nucleobases. After 5 ns fluctuations, compound 16 converged to a stable binding conformation and leveled off the rest of MD simulations time. This demonstrated that after 5 ns simulation and initial fluctuations in the all atom scale of ligand RMSD, the ligand acquired an equilibrium state characterized by the RMSD outline. At first glance, the complex between compound 18 and DNA minor groove seemed to be unstable. But a closer look at RMSD fluctuations revealed that the changes in RMSD were very small by the mean 0.066 ± 0.008 nm.
More informative illustration of the system stability could be achieved by computing the center of mass (COM) distance between ligand and DNA. To get the DNA-ligand interactions, we calculated the COM distance between ligand and DNA during the simulation time. The results are depicted in Figure 7.
Compound 16 showed dramatic changes in COM distance during simulation time. COM distance between the intercalated ligand (16) and DNA nucleobases was stable at mean distance 0.60 ± 0.05 till 4.54 ns when the first increase occurred. As can be seen in Figure 8, the new state seemed to be unstable and rapidly returned to the mean value. This trend repeated over the whole simulation time. To discover the system condition in these situations, complex structure between compound 16 and DNA was extracted 12 Soaking of DNA-ligand system was carried out by 100 ps constrained NVT followed by 100 ps NPT ensembles.
Following system equilibration, 50 ns MD production was performed in each case without any constraint.
Total energy, temperature and root mean square deviation (RMSD) were assessed to confirm the stability of trajectories. The average temperature during MD simulations at 300 K was found to be 299.0 K (± 3.6) and 298.2 K (± 2.3) in the case of models 1D32 and 102D, respectively. The temperatures of the both systems were stable for the whole simulation time. Fluctuation of temperature and energy is depicted in Figure 5. Obtained results showed that the conservation of energy over 50 ns MD simulations was convinced in both systems.    As can be seen in Figure 7, both ligands had overall stability in binding state. Compound 16 bound to DNA through intercalation between DC and DG nucleobases. After 5 ns fluctuations, compound 16 converged to a stable binding conformation and leveled off the rest of MD simulations time. This demonstrated that after 5 ns simulation and initial fluctuations in the all atom scale of ligand RMSD, the ligand acquired an equilibrium state characterized by the RMSD outline. At first glance, the complex between compound 18 and DNA minor groove seemed to be unstable. But a closer look at RMSD fluctuations revealed that the changes in RMSD were very small by the mean 0.066 ± 0.008 nm.
In Silico Study of Indenoquinoline diones from trajectories for each 5 ns (Figure 8).
We postulated that compound 16 was probably capable of separating DNA strands from each other. The rise in COM distance might be related to the departure of DNA strands from initial position. This movement increased the COM distance between compound 16 and DNA with regard to their initial conformation ( Figure  8A). The separation of DNA strands took place at the end of strand and we could not predict the Compound 16 showed dramatic changes in COM distance during simulation time. COM distance between the intercalated ligand (16) and DNA nucleobases was stable at mean distance 0.60 ± 0.05 till 4.54 ns when the first increase occurred. As can be seen in Figure 8, the new state seemed to be unstable and rapidly returned to the mean value. This trend repeated over the whole simulation time. To discover the system condition in these situations, complex structure between compound 16 and DNA was extracted from trajectories for each 5 ns (Figure 8).  We postulated that compound 16 was probably capable of separating DNA strands from each other. The rise in COM distance might be related to the departure of DNA strands from initial position. This movement increased the COM distance between compound 16 and DNA with regard to their initial conformation ( Figure 8A). The separation of DNA strands took place at the end of strand and we could not predict the situation in which the compound intercalated in the middle of DNA strands. In the case of compound 18, three increases followed by relevant decreases in COM distance could be detected. The fluctuations occurred between 15 and 22 ns, 27 and 35 ns and finally at 47 ns and continued to the end of simulation. Double helix structure was stable during simulation so this time the movement of ligand could be responsible for observed changes. Compound 18 moved along the minor groove. In this case an interesting fact was the returning of ligand into its initial position. Analyzing the binding energy between compound 18 and DNA revealed that initial complex had maximum binding energy ( Figure 9B). VDW interaction energy as a main contributor to total energy decreased due to the movement of ligand from its initial position. This result demonstrated that compound 18 preferred binding to DT reach domain. For more clarification, the contribution of each nucleobase in total binding energy is illustrated in Figure 10.
As it was expected, compound 18 attached to the minor groove of DNA mainly through VDW interactions. Moreover; compound 18 could favorably bind to DT rich domain through free movement in the minor groove.

Conclusion
The principal objective of the present study was to propose a binding model for cytotoxic indeno [1, 2-b] quinoline-9,11-diones as DNA VDW interaction energy as a main contributor to total energy decreased due to the movement of ligand from its initial position. This result demonstrated that compound 18 preferred binding to DT reach domain. For more clarification, the contribution of each nucleobase in total binding energy is illustrated in Figure 10. As it was expected, compound 18 attached to the minor groove of DNA mainly through VDW interactions. Moreover; compound 18 could favorably bind to DT rich domain through free movement in the minor groove.  VDW interaction energy as a main contributor to total energy decreased due to the movement of ligand from its initial position. This result demonstrated that compound 18 preferred binding to DT reach domain. For more clarification, the contribution of each nucleobase in total binding energy is illustrated in Figure 10. As it was expected, compound 18 attached to the minor groove of DNA mainly through VDW interactions. Moreover; compound 18 could favorably bind to DT rich domain through free movement in the minor groove.   interfering agents. Docking results indicated that these compounds could reversibly attach to DNA oligonucleotides with varied free binding energies ranging from about -6 to nearly 10 kcal.mol -1 . Docking studies confirmed that the incorporation of electron-withdrawing substituents such as chlorine atoms to the imidazole moiety of indeno[1,2-b]quinoline-9,11-dions might be associated with higher cytotoxic effects MD simulation studies confirmed the stability of docked complexes and showed that probable binding mode of indeno [1, 2-b] quinoline-9,11-diones within DNA double helix could be intercalation between DC3, DC5, DG4 and DG6 nucleobases and indene, imidazole, indene and quinoline moieties of the ligands, respectively (PDB ID: 1D32). Linear regression analysis showed that cytotoxic effects of aromatized scaffolds in Jurkat cell lines may not be related to the DNA intercalation. MD simulations also demonstrated that compounds bearing two methyl substituents at their quinoline ring could attach to the minor groove of DNA by VDW interactions with DT rich domain through free movement in the minor groove. N-phenyl substituted indeno [1, 2-b] quinoline-9,11-dions contributed to additional hydrophobic interactions with DC7 and DG8 (PDB ID: 1D32) and DA18, DT8, DT9 and DT19 (102D) nucleobases of DNA. The applied method can be extended to other biological systems involving DNA-ligand interactions. (4)